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The gas-liquid phase transition of the three-dimensional Lennard-Jones particles system is stud- 
ied by molecular dynamics simulations. The gas and liquid densities in the coexisting state are 
determined with high accuracy. The critical point is determined by the block density analysis of 
the Binder parameter with the aid of the law of rectilinear diameter. From the critical behavior 
of the gas-liquid coexsisting density, the critical exponent of the order parameter is estimated to 
be /3 = 0.3285(7). Surface tension is estimated from interface broadening behavior due to capillary 
waves. From the critical behavior of the surface tension, the critical exponent of the correlation 
length is estimated tohe v = 0.63(4). The obtained values of /3 and v are consistent with those of 
the Ising universality class. 



I. INTRODUCTION 

University and scaling are key concepts in the study 
of critical systems [UH , including liquid-gas systems [j- 
[1] , Ising model [I,l3l7 percolation model [l-|l2j , dimer 
model [l3|, etc. According to modern theory of critical 
phenomena [11, Q , critical systems can be classified into 
different universality classes such that the systems in the 
same class have the same set of critical exponents. 

The Lennard- Jones (LJ) particle system has been ex- 
tensively studied since it is the basic model of the off- 
lattice systems with short-range interactions. It is widely 
believed that the gas-liquid phase transition of the parti- 
cle system with short-range interaction belongs to Ising 
universality class IJ] and the LJ-particle system is no 
exception. However, there is no definite proof that the 
universality class of the liquid-gas phase transition of the 
LJ system is same as that of the Ising universality. In 
order to study the universality class of the LJ system, the 
precise information of the phase diagram, especially the 
location of the critical point, is required. Aside from the 
critical behavior, the determination of the phase bound- 
ary itself is important to study nonequilibrium phenom- 
ena such as bubble nuclei, etc [l^ . 

Recently, due to the advance in the computational 
power, many researchers have tried to determine the 
value of the critical exponents as accurate as possible, 
in order to confirm whether the universality class of the 
LJ system is really the same as that of the Ising model. 
Panagiotopoulos proposed the Gibbs ensemble Monte 
Carlo (GEMC) method to calculate the phase equilibria 
of the coexisting phase [l^ . The GEMC method involves 
two simulation boxes, one is in liquid phase and the other 
is in gas phase under a given temperature. There arc 



three types of Monte Carlo (MC) steps; displacements of 
particles, changing volumes of simulation boxes, and ex- 
changing particles between the simulation boxes. When 
the two systems reach their equilibrium states, wc can 
obtain the physical quantities such as densities of gas 
and liquid phases, vapor pressure and chemical potential 
of the coexisting state in the given temperature. While 
this method is simple, the system can be unstable near 
the criticality, e.g., the simulation box may change its 
phase from gas to liquid and vice versa because of the 
fluctuation. 

MoUer and Fischer proposed a different approach, NpT 
and test particle method [l7[. This method involves 



NpT simulation, i.e., the isobaric and isothermal condi- 
tion. The chemical potential is observed as a function of 
temperature by Widom's test-particle method [l^, and 
the critical point is determined by the crossing point 
of the chemical potential of liquid and gas branches. 
This method is more stable than GEMC near the crit- 
ical point; since only the volume fluctuates while both 
density and volume fluctuate in GEMC. Additionally, 
Okumura and Yonezawa proposed NVT and test 
particle method which is more stable than the last two 
methods since this method does not involve change in vol- 
ume. However, the calculation of the chemical potential 
with Widom's test-particle method is quite expensive, 
especially for larger systems. Wilding and Bruce com- 
bined the histogram reweighting (HR) techniques with 
the finite-size scaling (FSS) analysis [20|. The main idea 
of this method is to locate the apparent critical point 
from the matching condition between the ordering oper- 
ator distribution and that of Ising model with the help 
of the mixed-field theory. Perez- Pellitero et al. compared 
the HR techniques with the direct estimation of the crit- 



ical point using Binder parameter [21 
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While the HR techniques combined with FSS analysis 
are efficient to study the coexisting phase, the distribu- 
tion of the order parameter of Ising model at the critical 
point is required for the reference. Additionally, most of 
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past studies determined the critical point with the help 
of the knowledge of the Ising universality, and therefore, 
it is difficult to determine whether the gas-liquid phase 
transition of LJ-particle system belongs to the Ising uni- 
versality class. 

In the present study, we propose a simple method to 
study the gas-liquid coexisting phase using the molecu- 
lar dynamics (MD) simulations. The main purpose of 
our method is to determine the critical properties of the 
LJ gas-liquid system without a priori knowledge of the 
Ising universality, such as the values of the critical expo- 
nents, the value of Binder parameter at the critical point, 
or the distribution of the order parameter, etc. The crit- 
ical exponents of the order parameter is estimated from 
the liquid-gas coexisting density and the that of the cor- 
relation length is estimated from the interface thickness 
taking into account the effect of the capillary waves. We 
have found that the critical exponents of the order pa- 
rameter and correlation length are /3 = 0.3285(7) and 
V = 0.63(4), respectively, which are consistent with those 
of the Ising universality class [l|, [23 - |24j . 

This paper is organized as follows. In Sec. II, we de- 
scribe our model system and details of the simulation 
procedure. We first describe how to determine the coex- 
isting curve in Sec. II A, then we describe how to locate 
the critical point using the Binder parameter with the aid 
of the law of rectilinear diameter in Sec. II B. In Sec. Ill, 
we present results of our simulation data and critical ex- 
ponents determined from such data. Sec. IV contains our 
conclusions and discusses problems for further studies. 



-3.7 



P -3.8 



o 

D. 



-3.9 





Thermalization 


Observation 




Potential Energy 




Temperature 









2000 



4000 6000 
Time 



-I 0.8 

0.6 a 

2 

Q. 

0.4 1 
I- 

0.2 




8000 10000 



FIG. 1: Relaxation behaviors of potential energy (solid line) 
and temperature (dashed line). The system size is 64x64x 128 
with the number of particles 180288. Aimed temperature is 
0.9. While the relaxation of temperature is quick, that of the 
potential energy is slow. After thermalization, the heatbath 
is turned off. 



II. MODEL SYSTEM AND SIMULATION 
PROCEDURE 

A. Model system and coexisting curve 

We use the truncated Lcnnard-Joncs potential of the 
form 



V{r) = 4e 



C2 



Co 



(1) 



with the well depth e and atomic diameter a [25|. The 
coefScients C2 and cq are determined so that Vijc) = 
y'{fc) = with the cut-off length rc , i.e., the values 
of potential and force become continuously zero at the 
truncation point and V{r) =0 for r > rc- The two addi- 
tional terms arc introduced for good conservation of the 
total energy. The extra terms, which makes the force 
continuous at the truncation point, are required to per- 
form the MD simulations stably for long time and it is 
unnecessary for MC simulations. In the following, we use 
the physical quantities reduced by a, e, and the Boltz- 
mann constant fee, e.g., the length is measured with the 
unit of cr, and so forth. We set the cutoff length as 
Tc — 3.0. The simulation box is a rectangular paral- 
lelepiped with sizes y~ Ly x L^- The periodic bound- 
ary condition is taken for all directions. We set the ratio 
oi Lx '. Ly : Lz ~ \ '■ 1 '■ 2 ii\ order to achieve the 
gas-liquid coexisting phase. The studied system sizes are 
Lx = 32, 64, 96, and 128. We setup the initial configura- 
tion by placing particles densely in the region z < Lz/2 
and sparsely in the region z > Lz/2, and we equibrate 
the system with the heatbath. Nose-Hoover heatbath is 
used for controlling temperature [26| . Time integration 
scheme is the second order RESPA [23| for isothermal 
simulations and the second order symplectic integrator 
for isoenergetic simulations. The time step is chosen to 
be 0.005 throughout the simulations. We check whether 
the system has reached the equilibrium state by observ- 
ing the potential energy, since the relaxation of the po- 
tential energy is much slower than that of the tempera- 
ture (see Fig. [T]) . The potential energy U {t) behaves as 
U{t) — Uoq oc exp(— t/r) with the equilibrium value [/cq 
and the relaxation time r. We perform thermalization at 
least ten times longer than the relaxation time. Then we 
turn off the heatbath and start observation of the density 
profile at and near gas-liquid phase interface and surface 
tension. The simulations are performed with MPI paral- 
lelized MD program [2^ [2^ . The largest run, involving 
1583504 particles and 5000000 steps, takes about 8 hours 
using 1024 cores at SGI Ahix ICE 8400 EX. 

The functional form of the density profile at and near 
the gas-liquid interface is of the form, 

p{z)=^-^^^{t^nh[{z,^z)/X] + l) + p„ (2) 

where pi, pg, and A are the density of liquid, the density 
of gas, and the thickness of the interface. The typical 
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FIG. 2: Density profile at and near tiie gas-liquid interface. 
The conditions of the simulation are same as Fig. [T] The open 
circles are the simulation results and the solid line denotes 
the fitting result via Eq. ((2]) with parameters pi = 0.6729(3), 
Pg = 0.0397(3), Zc = 62.161(4), and A = 1.956(7). Statistical 
errors are smaller than the symbol size. 

density profile is shown in Fig. [2] By fitting Eq. ([2]) 
to the simulation results, we obtain the gas density pg, 
the liquid density pi, and the thickness A as functions of 
temperature. The thickness A is associated with the cor- 
relation length, and becomes longer as the temperature 
approaches to the critical point. The finite-size effect can 
be ignored while A is sufficiently shorter than the system 
size Lx- This means that larger system size allows us to 
approach closer to the critical point. Typical configura- 
tions of particles below, near, and above the critical point 
are shown in Figs. |3Ia), |31[b), and|3Uc), respectively. 



B. Surface tension 

When the phase interfaces are normal to the .z-axis, 
the surface tension 7 is defined by the following equation 

2,LxL, . |dr (...(r) - ^^:.iW±I-W) , (3) 

where tTxxj T^yy and n^z are the diagonal components of 
the stress tensor [sO]. The integral taken over the whole 
volume of the simulation box. The factor 2 of the l.h.s. 
comes from the fact that there are two interfaces in the 
system due to the periodic condition. For the particles 
interacting via two body potential V{r), Eq. ([3|) reduces 
as 

2^LxLy=j:^MlAzlfkv'^^^l (4) 

where is the distance between particles i and j , V'{rij) 
is the derivative of V{rij) with respect to r^j, Xij is the x 



component of the relative position vector between parti- 
cles i and j, and so forth. The sum is taken over all pairs. 
Note that, the kinetic term of the stress tensor is omitted 
since the kinetic term does not contribute to the surface 
tension. After the thermalization, the surface tension is 
sampled every 1000 steps and averaged. The observation 
is performed for 100000 steps for each temperature. 

C. Critical behavior 

The order parameter of the gas-liquid transition is the 
density difference Ap = p\ — pg. The order parameter 
becomes zero at the critical temperature Tc, and the be- 
havior near the criticality is 

Ap^e^ (5) 

with the reduced temperature e = \Tc — T\/Tc and the 
critical exponent j3. The surface tension 7 and the thick- 
ness of the interface A also show critical behavior as fol- 
lows, 

7 - (6) 

A ^ £-^ (7) 

with the critical exponent of the correlation length u. It 
is difficult to determine the critical temperature and the 
critical exponents simultaneously only from the above 
scaling relations without assuming the Ising universal- 
ity. Therefore, we determine the critical point using the 
Binder parameter [sij . 

In order to compute the fluctuation of density from 
NVT or NVE simulations, we adopt the block distribu- 
tion analysis [s^ . The main idea of this method is to 
divide the system into small cells, and each cell exhibits 
approximately the grandcanonical ensemble, while cells 
are not completely independent from each other. Sup- 
pose there are a cubic system with linear dimension L 
containing N particles. The simulation box with linear 
dimension L is divided in to cells of size L/il/f,, where 
Mb is an integer. Each cell contains different number of 
particles, and therefore, has different density pi, i being 
the label of the cells. The moments of density required 
for calculation of Binder parameter are defined to be. 




= j^T.(p^-p)'^ (9) 
- -^T.(p^-p)'- (10) 

^ i 

Using the above definition, we can define the Binder pa- 
rameters as follows: 

U{M,) = (11) 
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FIG. 3: Typical configurations of particles. Small systems 
with sizes {Lx,Ly,Lz) = (32, 32, 64) are shown for the visibil- 
ity, (a) Below critical points (T = 1.00) with 18432 particles. 
The gas (left) and liquid (right) phases are clearly separated, 
(b) Near critical points (T = 1.09) with 18824 particles. The 
correlation length becomes longer and the phase interface be- 
comes unclear, (c) Above critical points (T = 1.11) with 
17576 particles. The system becomes uniform. 



where (• • ■ ) denotes ensemble average. With different di- 
viding number Mf,, we can obtain the Binder parameters 
for different volume sizes from a single run for a suffi- 
ciently large system. 

The crossing point of the Binder parameter gives the 
critical point. Since we perform NVT ensemble, the den- 



sity and the temperature should be given simultaneously. 
Therefore, we assume the law of rectilinear diameters [33| 
that the average density shows linear dependence on tem- 
perature as (pg + p\)/2 = aT + b. While some fluids 
exhibit strong deviations from the above linear relation, 
simple pure fluids such as oxygen and xeon follow this 
law within high accuracy [s^ . The studied system is a 
cube with linear dimension L = 128. First the system 
is thermalized for 200000 steps with heatbath, and the 
Binder parameters are observed every 1000 steps. The 
observation is performed for 100000 steps. The Binder 
parameters are calculated for Mf, = 16, 24, and 32. 

While the critical exponent can be also estimated from 
the finite-size analysis of Binder parameter, it is better to 
estimate the critical exponent /3 from the liquid-gas co- 
existing density. Because the coexisting density can be 
determined more accurate than Binder parameter, since 
the coexisting density is the first order moment of the 
density while Binder parameter contains the higher mo- 
ments. Therefore, we use Binder parameter only for lo- 
cating the critical point and we estimate the critical ex- 
ponent from the liquid-gas coexisting density. 

We first calculate the phase diagram, and determine 
the coefficients of the linear relation. Then we calculate 
the Binder parameter on the line in order to determine 
the critical temperature. Using the obtained critical tem- 
perature, we determine the critical exponents. 
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FIG. 4: (Color online) Phase boundary between the gas 
and the liquid phase obtained by the simulation with Lx ~ 
32, 64, 96, and 128. Statistical errors are smaller than the 
symbol size. The experimental result of neon is denoted by 
"Exp" [35I ]. The finite-size effect is hardly observed. The 
dashed line denotes the law of rectilinear diameters and the 
solid circle denotes the critical point obtained from the anal- 
yses of the Binder parameter. 
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FIG. 5: (Color online) Temperature dependence of the 
Binder parameter with different block sizes. The system is di- 
vided in to cells of size L/Mi,, where L = 128 and Mi, — 16, 24, 
and 32. As the value of Mb becomes larger, the size of cells 
becomes smaller. From the intersection point, the critical 
temperature is estimated to be Tc — 1.100(5). 
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FIG. 6: (Color online) Order parameter Ap = pi— pg is shown 
as a function of the reduced temperature e = \Tc — T\/Tc. The 
decimal logarithm is taken for the both axes. Statistical errors 
are smaller than the symbol size. Significant differences are 
hardly observed between the data of different system sizes. 
The critical exponent are determined to be /? = 0.3285(7). 
The solid line is the eye guide with the slope /3 — 0.3285. 

III. RESULTS 

The gas-liquid phase boundary obtained from the sim- 
ulations are shown in Fig. 21 The obtained gas and liquid 
coexisting densities do not exhibit finite-size effect, i.e., 
the values of different system sizes are within the sta- 
tistical errors. We also include the experimental data 
of neon which is denoted by "Exp" [S^- The obtained 
results for the largest systems — 128 are listed in Ta- 
ble, m We confirm that the law of rectilinear diameters 
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TABLE I: The data obtained from the simulations with the 
size 128 X 128 x 256. The listed quantities are temperature 
T, the number of particles A'', the gas density pg, the liquid 
density pi, the surface tension 7, and the interface thickness 
A, respectively. 



(pg + pi)/2 = aT + b is well satisfied. We determine 
the coefficients to be a = —0.195(1) and b = 0.533(1). 
The obtained line is shown as the dashed line in Fig. S) 
Note that, we extend the line into the supercritical phase 
while the law of rectilinear diameters does not make sense 
outside the coexisting phase, since we want to calculate 
binder parameter in the region over the critical point in 
order to have the intersection point. 

The Binder parameters are calculated on this line and 
they are shown in Fig. [5] The Binder parameter does not 
intersect exactly at one point as reported before [2lj , but 
we can determine the critical temperature with reason- 
able accuracy as Tc = 1.100(5). Using the critical point, 
we determine the value of the critical exponent /? from 
the slope of In A/o/ Ine as shown in Fig. |6l The obtained 
value (3 = 0.3285(7) is consistent with the reported value 
/3 = 0.3269(6) of the spontaneous magnetization [l[ of the 
three-dimensional Ising model on the simple cubic lattice 
obtained byTalapov and Blote HI or /3 = 0.325(5) by 
Ito, et al. [23] of the Ising universality class. The critical 
behavior of the surface tension is shown in Fig. [71 The 
critical exponent of the correlation length is determined 
to he V = 0.640(4) which is slightly larger than that of 
the Ising universality class v = 0.63002(10) 

We also study the critical behavior of the interface 
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FIG. 7: (Color online) Temperature dependence of the surface 
tension. The decimal logarithm is taken for both axes. The 
critical exponent is determined to be = 0.640(4). The solid 
line is the eye guide with the slope 2u = 1.28. 



thickness. We find that the interface thickness exhibits 
strong finite-size effect as shown in Fig. |8] (a), while the 
finite-size effect on the the surface tension is almost neg- 
Ugible. These finite-size dependencies come from the in- 
terface broadening due to the capillary waves [H, [STf . 
According to the capillary wave theory, the system size 
dependence of the interface thickness A is expressed to 
be, 



C 



T 

47 



InL, 



(12) 



where C is a system-size independent constant for a fixed 
temperature. Note that, the temperature T and system- 
size L used here arc dimensionless values. The system- 
size dependence of the interface thickness for T = 0.9 
is shown in the inset of Fig. [S] (a). From this finite-size 
dependence, we estimate the surface tension as a func- 
tion of temperature. The critical behavior of the surface 
tension obtained from the interface thickness is shown in 
Fig. [H] (b) . The critical exponent of the correlation length 
is determined tohe v = 0.63(4) which is consistent with 
the value of the Ising universality class. 



IV. DISCUSSIONS 

Wc determine the phase boundary between gas and 
liquid phases by achieving gas-liquid coexisting state 
with MD. This method is simple and provides the phase 
boundary between liquid and gas phases with high accu- 
racy. Since we have gas-liquid interface in the system, we 
can check the finite-size effect on the obtained results; if 
the linear length of the system is much longer than 
the interface thickness A, then the coexisting densities of 
liquid and gas are expected to be free from the finite- 
size effect. In our experience, the finite-size effect on 
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FIG. 8: (Color online) (a) Temperature dependence of the 
surface thickness. The significant finite-size effect is observed. 
(Inset) The interface broadening for T = 0.9. Using Eq. ([T2l) . 
we can estimate the surface tension, (b) Temperature de- 
pendence of the surface tension obtained from the interface 
broadening behavior. The critical exponent is determined to 
he V = 0.63(4). The solid line is the eye guide with the slope 
2u = 1.26. 



the coexisting density is vanished when Lz > lOA. We 
determine the critical temperature from the block den- 
sity analysis of the Binder parameter, which allows us to 
perform finite-size analysis only from one simulation run 
with large system size. Additionally, this method does 
not require the distribution of order parameter of Ising 
model at the criticality, and therefore, we do not have to 
consider the mixed-field theory. 

The critical exponents /3 and are determined from the 
behaviors of the order parameter and the surface tension. 
While the critical exponents of the order parameter /S is 
consistent with the Ising universality class, the exponent 
of correlation length v obtained from the surface tension 
is found to be slightly larger than that of the Ising univer- 
sality. While this deviation can be due to the capillary 
waves, it is difficult to remove the effect since the sur- 
face tension does not exhibit apparent finite-size effect. 
Potoff and Panagiotopoulos studied the critical bchav- 
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ior of the surface tension using HR technique and FSS 
analyses, and reported v = 0.71(4) which is much larger 
than the Ising universality However, they did not 

consider the effect from the capillary waves. Considering 
the effect of the capillary waves, we obtained the criti- 
cal exponent of the correlation length v = 0.63(4) which 
is consistent with the Ising universality. Note that, the 
statistical error became larger than that of the value ob- 
tained directly from the stress tensor using Eq. (|4]) . This 
is because the surface tension analysis with considering 
the capillary waves requires two steps; first measuring 
the surface thicknesses for various system sizes and then 
extracting the surface tension from its size dependence. 
While the obtained values of exponents are consistent 
with that of the Ising universality class, it is still difficult 
to declare that the universality class of the LJ system 
belongs to the Ising universality since the accuracy is in- 
sufficient especially for the exponent of the correlation 
length. Some techniques are required which can esti- 
mate the surface tension accurately. Also the distances 
from the criticality is insufficient enough to discuss the 
universality class precisely. In order to get closer to the 
critical point, larger and longer simulations are necessary. 
It should be one of the further issues. 

Some of the past studies reported the effect from the 
correction to scaling [l^, [2l| . They first obtained the ap- 
parent critical temperatures for several system size, and 
they took the infinite size limit considering the correc- 
tion. However, we did not find any sign of the correction 
to scaling in Fig. [HI It is difficult to determine the value of 
the correction to scaling exponent if it exists. While the 
obtained values of critical exponents are consistent with 
the Ising universality without considering the correction. 



it would become more important in order to investigate 
more accurately. 

Very recently. Ma and Hu [H-Iil] used MD to study 
relaxation and aggregation of polymer chains, in which 
neighboring monomers of a polymer chain are connected 
by rigid bonds [1^ |4l| or springs with various spring con- 
stants [13, 112] . They found that when the bending-angle 
depending potential and the torsion-angle depending po- 
tential are zero or very small, polymer chains tend to 
aggregate [4l| - |43j . Such results are useful for understand- 
ing mechanism of protein aggregation [4l|, |42| related 
to neurodegenerative diseases, including Alzheimer's dis- 
ease (AD), Huntington's disease (HD), Parkinson's dis- 
ease (PD), etc. One may argue that polymer chains can 
aggregate only when they are inside the phase bound- 
ary. It is of interest to extend the method of the present 
paper to calculate the phase diagram to calculate phase 
diagrams of polymer chains of various chain lengths [3l ■ 
Such phase diagrams will be useful for understanding un- 
der what temperature and density conditions, polymer 
chains would aggregate. 
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